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INVERSE PROBLEMS IN RADIATIVE 
TRANSFER: A REVIEW 


by 

Barney J. Conrath 


ABSTRACT 

An elementary introduction to the inverse problem of radiative 
transfer as applied to remote measurements of infrared radiation 
from planetary atmospheres is presented, along with a review of 
work which has been done on the problem. Particular attention is 
given to the problems of inferring vertical temperature profiles 
and water vapor distributions in the earth's atmosphere. The 
principle methods which have been de veloped for solving these 
problems are discussed briefly. Examples of applications of tem- 
perature and water vapor inversions to both synthetic data and data 
taken with an IRIS instrument at balloon altitudes are presented. 
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INVERSE PROBLEMS IN RADIATIVE 


TRANSFER: A REVIEW 


INTRODUCTION 

General 


In the direct problem of radiative transfer as applied to a planetary at- 
mosphere, the vertical structure of the atmosphere is specified along with the 
energy sources, and the radiation field is calculated. In particular the spectral 
intensity of the upwelling electromagnetic radiation at the effective top of the 
atmosphere can be predicted. In the inverse problem the outgoing intensity is 
assumed to be obtained observationally, and an effort is made to infer informa- 
tion on the vertical atmospheric structure. In the present paper we shall con- 
sider a narrowly restricted class of the inverse problem; inference of vertical 
temperature profiles and vertical constituent distributions from remote meas- 
urements of thermally emitted radiation. 

The primary purposes of the paper are to provide an introduction to the 
subject and summarize the present level of development. The problem of in- 
ferring the vertical temperature profile will first be considered, and a review 
of the basic methods which have been developed will be given. The important 
questions of the accuracy of the inferences will be discussed, and examples of 
applications to both synthetic and real data will be presented. Finally we shall 
discuss some recent work on the problem of remotely inferring the distribution 
of water vapor in the lower terrestrial atmosphere. 

Physical Description of the Problem 


As an example, let us consider the problem of inferring the temperature 
profile in the lower atmosphere of the earth by observing the spectral intensity 
in the 15m CO 2 absorption band with an instrument located above the effective 
upper boundary of the atmosphere. The observations consist of a spectral scan 
from the band center out into the band wings. From a knowledge of the ab- 
sorbtivity of C0 2 as a function of wavenumber, we can obtain the relative 
optical depths. With the additional information that C0 2 is uniformly mixed 
with a known mixing ratio in the part of the atmosphere we are considering, 
we can relate the optical depths at the various wavenumbers to actual pressure 
levels in the atmosphere. In other words we can calculate the transmissivity 
of the atmosphere from a given pressure level to the top as a function of 
wavenumber. 
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In Figure 1 we have plotted the pressure level at which the transmissivity 
drops to 1/e. This figure corresponds to a resolution equivalent to 0.1 cm -1 
which is beyond the reach of present satellite instrumentation, but serves to 
illustrate our point. In scanning from the band center toward the wing we see 
progressively deeper into the atmosphere. Thus, from measurements at se- 
lected wavenumbers we should expect to be able to infer information on the at- 
mospheric temperature at a number of different pressure levels, and in this 
way piece together a temperature profile. 

Now the intensity at a given wavenumber, instead of depending on the tem- 
perature at a single level, actually contains contributions over* a number of 
levels on either side of the e-folding pressure. To isolate the temperature at a 
given level we would have to subtract out contributions from all other layers 
which is equivalent to differentiation of the data. Thus, we might intuitively 
expect to encounter problems with the stability of our solutions. 

Mathematical Formulation 


The solution to the equation of radiative transfer for a nonscattering 
planetary atmosphere in local thermodynamic equilibrium can be written in 
the well known integral form 


1(10 - I,(i/)r(v, X s ) - B[v, T(x)] 


f 

•'o 


d'T(v, 

dx 


x) 


dx 


(1) 


where I(^) is the spectral specific intensity emerging at the top of the atmosphere, 
B(^, T) is the Planck intensity at wavenumber v and temperature T, x is any 
independent variable which is a single valued function of pressure, and t(v, x) is 
the transmissivity at wavenumber v of the column of atmospheric gas between 
level x and the top of the atmosphere. The subscript s refers to the planetary 
surface. For simplicity, we shall assume that our measurements are from 
sufficiently opaque parts of the spectrum that the boundary term in (1) can be 
neglected. 

We see from (1) that the outgoing intensity can be interpreted as a weighted 
mean of the Planck functions associated with each layer of the atmosphere. The 
weighting function is just dr/dx. Typical weighting functions for the 15 ^ C0 2 
band are shown in Figure 2. These functions were calculated by Kunde (private 
communication) for an instrumental response function 5 cm -1 in width. From 
this figure we can obtain an idea of the relative contribution of each atmospheric 
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layer to the outgoing intensity at each wavenumber. The "width" of the weighting 
functions is determined primarily by the way in which the optically active gas is 
distributed and for pure monochromatic radiation is ~1.4 atmospheric scale 
heights at the 1/ e points. Averaging over a finite spectral interval in which the 
absorption coefficient is changing will tend to increase this width somewhat. 

The spectral resolution required depends on the number of pieces of informa- 
tion we are trying to retrieve and the levels of the atmospheres at which in- 
formation is sought. In general by increasing the resolution we can reach 
higher levels in the atmosphere. A resolution equivalent to about 5 cm -1 seems 
to be the lower practical limit for inferring five or six pieces of information in 
the troposphere and stratosphere. Relatively broad band measurements could be 
used to infer limited information on the troposphere. 

The temperature inversion problem can now be stated as follows; given 
I(iv) anddr/dx, solve (1) for T(x). In its present form, (1) is a nonlinear integral 
equation in T(x). The problem can be simplified however by making use of the 
fact that the Planck function is a slowly varying function of wavenumber in the 
spectral region and at the temperatures we are considering. Hence we can 
linearize (1) by using an approximation of the form 


B[v,T] * a(v) + b(*)B[v 0 , T] 


( 2 ) 


where is some reference wavenumber near the center of our spectral meas- 
urements. Letting B [^ 0 , T(x)] 3 B(x), (1) becomes 


eO) - 



B( x ) K( v , x ) dx 


( 3 ) 


where 


g(^) 


I(V) - a(u) 
b(v) 


K(v, x) 


dr(v, x) 
dx 
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We can now regard B(x) as our unknown function since there is a one-to-one 
correspondence between B(x) and T(x). Equation (3) is a linear Fredholm in- 
tegral equation of the first kind. Such equations are notorious for their in- 
stability, the stability depending on the nature of the kernel K(v, x). 

We have already indicated the physical nature of the instability. Some in- 
sight into the mathematical nature of the instability can be gained from the 
following consideration. Suppose B 0 (x) represents the solution to (3). Now add 
to the solution a term of the form sin 2tt fx, and insert this into the right hand 
side of (3), carry out the integration and compare the result with g(i'). If the 
frequency f is sufficiently high so sin 277 fx changes rapidly compared to 
K(v, x) for all v considered, then adjacent negative and positive half cycles 
tend to cancel, and the contribution of this term will be small. Hence, the total 
integral will differ only slightly from g(^). Thus, we see that if g(v) is im- 
perfectly known as will be the case in practice, wildly oscillating solutions may 
result in place of the relatively smooth solution sought. 

The problem is further compounded by the fact that we will have available 
measurements corresponding to g(v) only at say N discrete values of v. Thus, 
even with perfect data, we cannot hope to obtain a unique solution without im- 
posing additional constraints upon the problem. In other words, we must specify 
an interpolation between the measured data points. 

We may summarize the problem as being one in which a solution B(x) to 
(3) is sought from an imperfect specification of g(^) at N discrete values of v. 
Any attempted solution must specify in effect an interpolation of g(^), and a 
means whereby the basically unstable solution can be objectively smoothed. 


METHODS DEVELOPED FOR TEMPERATURE INVERSION 
Linear Methods 


Perhaps the most direct way to attack the problem is to attempt to expand 
the Planck function B(x) in terms of some function set Fj (x). Such an expansion 
must of course be truncated at a number of terms equal to or less than the 
number of independent measurements we have available to us. On substitution 
of the expansion into (3) we obtain 


M /.x 

iO) - ^ a j J F. (x)K(i', x)dx 
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where the a^'s are the expansion coefficients to be determined. Thus, in 
choosing the set Fj (x) we have implicitly specified the functional form of the 
interpolation for g(^). If we like we can regard the process of finding the a - } 's 
as a curve fitting of (4) through the measured points g i = g(j'i). We obtain the 
set of linear equations 


Bi 



i - 1, 2, • • • N 
N > M 


( 5 ) 


where we have let 


A. . 


* j 



Fj (x)K(v i 


x) dx 


For the sake of simplicity we shall assume N = M, but what follows applies to 
the more general case also. It is convenient to write (5) in matrix notation 


g - A a 


( 6 ) 


We can now obtain a direct solution by solving (5) or (6) in the form 


a - A 1 g 


( 7 ) 


Such a solution applied to noisy data usually leads to catastrophic results, be- 
cause the matrix A is generally ill conditioned. This is equivalent to saying that 
the equations (5) do not possess a high degree of independence. This in turn 
can be traced back to the fact that the relative contributions from the various 
levels in the atmosphere are highly overlapping from one measured intensity 
to another. 

\ 
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The propagation of the random error in the measurements of g can easily 
be analyzed and can be expressed in the form 


Sa 

a 


'Xj 



( 8 ) 


where the error amplification factor G is in general dependent on the function 
set employed and the number of terms in the expansion. A typical value for a 
5 term expansion in terms of line segments is G ~ 10 4 . Obviously such an 
approach is impractical because of the extremely high demands placed on the 
accuracy with which g must be determined. In order to obtain a with an ac- 
curacy of only 10%, it would be necessary to reduce the experimental random 
error to less than one part in 10 s . 

Since a direct inversion of data containing what might normally be con- 
sidered reasonable errors leads in general to wildly oscillating solutions, it 
becomes necessary to introduce some type of smoothing process. Perhaps the 
most promising method of smoothing is that utilized by Wark and Fleming 
(1966) which employs a method of solution initially suggested by Phillips (1962) 
and extended by Twomey (1963). 

The method can be viewed essentially as follows. Assume we have N 
measurements g t with a certain probable error associated with them. Now let 
us allow ourselves the freedom of modifying each g i by some additive amount 
e i . Of course as soon as we do this we have available to us once again an in- 
finite number of solutions, and we must impose some additional constraints. 

One obvious constraint to impose is that the e.'s should not exceed the probable 
error of our measurements. We still have an infinite number of solutions 
available to us, but they are now members of a rather restricted set. From 
this set of solutions we pick that solution which minimizes some quadratic 
form of the a ( 's 


<3 = ». 

» . j 


( 9 ) 


where the coefficients H - t ■ are determined by the type of smoothing constraint 
which we wish to impose. The specific quadratic form chosen depends on our 
a priori knowledge of the form of the solution. For example we may wish to 
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minimize the mean square curvature of the solution or the mean square deviation 
from the mean or perhaps the mean square deviation from some trial solution 
based on a climatological mean temperature profile. 

The smoothed solution a (s) can be expressed in the form 


a ( s > - (A A + >H) 1 Ag 


( 10 ) 


where A is the transpose of A, H is a matrix whose elements are the coefficients 
H . . , and 7 is a Lagrange multiplier. The value of 7 can in principle be related 
to the probable error in the data. In practice it is established through trial and 
error. In general increasing 7 increases the amount of smoothing. For 7=0, 
(10) reduces to the direct inversion (7) and a( s) = a. Thus, if we like we can 
view ( 10 ) as a means of introducing a specified type of smoothing with the de- 
gree of smoothing being conveniently controlled by the value of a single scalar 
parameter. Arithmetically, the effect of increasing 7 can be viewed as simply 
making the matrix to be inverted more diagonally dominant and hence better 
conditioned. This follows from the fact that for most reasonable types of 
smoothing H is strongly diagonally dominant. 

An example of this type of smoothed solution is given in Figure 3 where 
data from a balloon borne IRIS instrument obtained near Palestine, Texas, 

8 May 1966 have been employed. The function set used consists of four line 
segments requiring the determination of five coefficients. The random error 
of the measurements is uncertain, but is estimated at about 2%. The curve 
marked 7 = 0 is the unsmoothed solution. The curve for a value of 7 giving 
the best fit is shown along with a solution employing a relatively large value 
of 7, resulting in an oversmoothing. The constraint employed was minimiza- 
tion of the second differences of the amplitudes at the break points. 

The question of the true information content of the data is an important 
one. We would like to ascertain how many independent pieces of information 
we can expect to infer from data of a given accuracy. Assume that using N 
intensity measurements we have successfully determined N coefficients of an 
expansion for the temperature profile. We can now substitute the inferred 
profile back into the equation of transfer and use the resulting relation for in- 
terpolation or extrapolation. Now assume we have an additional measurement 
available to us. We can use our extrapolation formula to predict the value of 
the additional measurement, and if the predicted value agrees with the measured 
value to within the probable error of the measurement, then we cannot hope to 
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improve on our solution by adding an additional term to the expansion. Indeed, 
any attempt at inferring an additional expansion coefficient can be expected to 
increase the instability. 

While we are limited to the number of independent pieces of information 
which we can obtain, additional information can be used in a redundant sense 
to improve the effective signal -to-noise ratio. To realize a significant im- 
provement, one must have available a number of additional measurements since 
the effective signal-to-noise ratio is approximately proportional to the square 
root of the redundancy. 

The number of independent pieces of information which can be inferred 
from a given set of data depends on the degree of independence of the kernel 
functions and the signal-to-noise ratio. Twomey (1966) has provided an 
elegent mathematical method of analyzing this problem. It appears that at most 
about 5 independent pieces of information can be inferred from present state- 
of-the-art satellite borne instruments. 

Expansion Functions 


Since we have a limited number of expansion terms available to us, it is 
important to choose a function set which gives us the most rapid convergence 
possible. Most of the classical orthogonal function sets have been tried, and 
of these the best appears to be trigonometric functions. Straight line segments 
can be made to fit well, but in general there is no good a priori way of choosing 
the break points. 

In the case of temperature profiles in the earth’s atmosphere, we have 
available to us a great deal of information on the statistical behavior of the 
profile for a given season and latitude. A convenient means of incorporating 
such information into our inversion technique is provided by the use of empirical 
orthogonal functions. Application of this technique to the temperature inversion 
problem has been made by Alishouse, et al. (1966). 

Let us assume we can treat the profile B(x) as a stochastic function whose 
covariance is given by 


R(x, x'.) = ([B(x) - <B(x)>] • [B(x')- <B(x')>]) (11) 


where the angular brackets denote ensemble averaging. It can be shown 
(cf. Obukhov, 1960) that the eigenfunctions of R(x, x') taken as an integral 
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operator 


Jr(x, x')<£ n (x')dx' - k n ct> n (x) 


( 12 ) 


form the most natural set for expansion of members of the ensemble in the fol- 
lowing sense. Let b(x) = B(x) - (B(x)) be approximated by N terms of an ex- 
pansion in an orthogonal function set and let this approximation be denoted by 
b (N) (x). Then define a mean square error for the ensemble by 


a 


2 

N 


J <(Jb(x) -b (N) (x)] 2 ^ dx 


(13) 


Of all possible sets of orthogonal functions, ct n 2 is a minimu m for the set <£ n (x) 
obtained from (12) and ordered by decreasing eigenvalues . Thus we would 
expect the set <t- n (x) to provide a good set for use in our inversion problem, 
providing we can accurately estimate R(x, x'). 

To gain some insight into the physical interpretation of the <£ n ’s we make use 
of the relation 


R(x, x ' ) = 

n - 1 


(14) 


which is derivable from (12) and the orthogonality properties of the <ps. R(x, x') 
is the relative correlation of the derivation from the mean at levels x and x ' . 

The first term in (14) is the first approximation to R(x, x') and contains the 
overall gross features of the function. Thus (x) reflects the large scale 
statistical relationship of the derivations from the mean. The higher order 
functions are measures of the correlation of finer structure. 

In Figure 4 we show two functions, the first order eigenfunction and the 
fifth order eigenfunction for an ensemble of twenty-five temperature profiles 
narrowly restricted in geographic location and time. The two functions shown 
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serve to illustrate the way in which the structure increases as we go to higher 
eigenfunctions. By looking at 4 > x , we see that for this ensemble when we have a 
tropopause cooler than the mean, it is statistically related to both surface 
temperatures and stratospheric temperatures warmer than the mean. 

As a numerical experiment we have chosen one member of the ensemble 
and calculated synthetic "data" for five points in the 15/x C0 2 band. These data 
were then used in an inversion employing a five term expansion in the empirical 
orthogonal function set. The results are shown in Figure 5. The solution dis- 
plays an amazing amount of structure for a five term expansion. We next con- 
sidered a profile which was not a member of the ensemble. The resulting inversion 
is shown in Figure 6 where we see the agreement is not as good this time, 
particularly at the tropopause. While the deviation from the mean of the second 
sounding was no greater than that for the first sounding, the behavior of the 
second at the tropopause was somewhat atypical which apparently accounts for 
the poorer fit. 

In employing empirical orthogonal functions, we are restricting ourselves 
in such a way that we can generally not hope to recover unusual profiles. Some 
profiles will be represented better than others. In other words, if we are going 
to employ statistical methods we must expect statistical results. The degree of 
success of the method depends on how accurately we can estimate R(x, x'). This 
is a subject of current investigation. 

Nonlinear Method of King 


We shall now consider a somewhat different approach to the temperature inver- 
sion problem developed by King (1964). In this method it is assumed that (1) 
can be written in the form 


I(s) 


f 


B(x) e~ 


dx 


(15) 


for a channel sufficiently opaque that the boundary term may be neglected. 

The parameter s is a function of wavenumber. Form (15) assumes the trans- 
mission function can be represented by the analytical formr(v, p) = exp[-s(^) x(p)]. 
This form is valid for purely monochromatic radiation, but is only approximate 
for transmissivities averaged over finite band widths. It has been found, how- 
ever, that transmissivities in the 15M C0 2 band can be fit fairly well by this 
analytic form over a limited wavenumber range. 
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If (15) is assumed to be valid, the outgoing intensity is formally just the 
Laplace transform of the Planck function profile. Unless I(s) can be expressed 
in some analytic form whose Laplace transform is known, this helps us very 
little since numerical inversion of Laplace transforms is in general an unstable 
process (Bellman, et al., 1966) as we might expect since we have done nothing 
to change the physical nature of the problem. Let us turn the problem around 
and assume some physically reasonable form for B(x), calculate I(s) and at- 
tempt to fit the resulting expression to the data just as we did in the linear 
method. 

It is convenient to integrate (15) by parts 


I(s) - B(0) + e sx dx (16) 

. Jo 


where B(0) is the Planck intensity at the top of the atmosphere. Let us assume 
we can represent B(x) as a set of n isothermal slabs so we can write 


dB(x) 

dx 


Z>i s ( x_x i) 

j = l 


(17) 


where AB j is the change in B at the slab boundary x- . Instead of specifying the 
Xj's a priori as is done in the linear methods when line segments are used, they 
are treated as unknowns which are to be determined from the measurements 
along with the ABj's. Substitution of (17) into (16) gives 


I(s) 


B( 0) + 



(18) 


Thus, we are requiring that our data be fit by a series of exponentials with the 
free parameters B(0) , n values of bj 's, and n values of AB j 's. If we have meas- 
urements of I at 2n + 1 discrete values of s, then (18) becomes a set of 2n + 1 
nonlinear equations in 2n + 1 unknowns . 
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There exists an algorithm for solving such a set of equations (cf. Lanczos, 
1964). This algorithm requires values of I(s) at integral values of s. Physically 
this means that data must be sampled at integrally spaced values of the absorp- 
tion coefficient. (To apply the algorithm directly, it is necessary to specify 
B(0) by some other means. This point will be discussed further below.) 

In order that our solution be physically meaningful, the x ^s must be real 
and positive. If we get one or more x ^'s which do not meet this criterion, then 
we must conclude that there is no atmosphere containing n isothermal slabs 
which can produce the values of I(s) at the measured points. This may be due 
either to noisy data or simply that n isothermal slabs are not an adequate 
representation of the true thermal profile. 

For the thermal inversion problem being considered, it has been found 
empirically that the AB . associated with a negative or imaginary x ■ is fre- 
quently one or two orders of magnitude smaller than the remaining ABj's. 

When this occurs, the term in (18) involving the unacceptable x, can be dropped 
and the resulting n - 1 isothermal slabs taken as a smoothed solution. 

In practice, the atmospheric temperature profile can be approximated some- 
what better by ramps than slabs. The use of ramps can be implemented by in- 
tegrating (16) by parts to obtain 


sl(s) 


sB(O) + B' (0) + 


f d 2 B 

J 0 dx 2 


dx 


(19) 


We assume the vertical profile can be approximated by 


d 2 B 
dx 2 


2^AB.' 8(x-x.) 
j= i 


(20) 


where the AB 's are the discontinuous changes in slope. This results in 


sl(s) 


sB(0) + B' (0) + 



( 21 ) 
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We still require only 2n + 1 measurements to determine the 2n + 2 parameters 
since the additional relation required is the s = 0 case which is physically 
equivalent to the requirement that the temperature be finite at arbitrarily large 
optical depths. The techniques described for the isothermal slab model can now 
be applied. 

When the nonlinear method was first developed, the upper boundary conditions 
were specified by what amounted to a guess. However, better techniques have 
now been developed. For example a differentiation of (21) with respect to s elim- 
inates B' (0) 


I(s) + si' (s) 



d 2 B 
dx 2 


e~ sx dx 


( 22 ) 


A second differentiation with respect to s eliminates B(0) 


si" (s) + 21' (s) 



d 2 B 
dx 2 


e 


- SX 


dx 


(23) 


In principle we could apply the algorithm directly to (23). In practice it is 
obviously highly undesirable to try to obtain I ' (s) and especially I " (s) from 
noisy data. The method currently being employed by King involves use of 
(22) with 2n + 1 measurements in the 15m C0 2 band and an additional measure- 
ment in the atmospheric window. The value of B(0) is chosen by an iterative 
procedure to make the temperature profile agree at the surface with the surface 
temperature inferred from the window measurement. 

Because of the nonlinear nature of this method, error analysis is more dif- 
ficult than in the linear method. The best approach seems to be application of 
the method to synthetic data from a wide variety of model atmospheres. Ran- 
dom noise can be applied to these ’’data" to simulate instrumentation effects. 

A program of this type is currently being carried out by King. Preliminary re- 
sults indicate that with realistic satellite instruments one can hope to infer only 
two ramps consistently and in rare instances three ramps. It should be noted 
that a two ramp atmosphere corresponds to five independent pieces of informa- 
tion, which is essentially the same conclusion reached for the linear method. 
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We must consider whether or not a two or three ramp temperature profile 
provides useful information for the earth's atmosphere. Figure 7 shows a 
three ramp inversion of a tropical-type atmosphere, and the agreement is found 
to be fairly good. However, a three ramp inversion of a midlatitude profile does 
not give a particularly good fit in the stratosphere as shown in Figure 8. The 
reason for this is simply that a profile of this type cannot be well represented 
by three ramps. 

The nonlinear method has the advantage of objective smoothing of noisy data 
in the sense that when noise results in a spectrum which is physically incon- 
sistent with a temperature profile of n ramps a spectrum resulting from a 
smaller number of ramps is substituted. It has the disadvantage that the ex- 
ponential kernel required is not capable of giving a perfect fit of actual trans- 
missivities. It is not possible to incorporate a priori knowledge of the earth's 
atmosphere into the formulation. The method may be better suited for applica- 
tions to radiometric data from the atmospheres of other planets where little 
a priori information is available. An example of an application to synthetic data 
for a model Martian atmosphere is shown in Figure 9. 

Temperature Inversion in Cloudy Atmospheres 


In the discussions given so far we have implicitly assumed that we are 
dealing with a completely cloud-free atmosphere. The formulations given above 
are equally applicable for situations in which a homogeneous cloud layer exists 
throughout the field of view of the measuring instrument. Cases will frequently 
occur in which the field of view is only partially cloud filled or contains cloud 
layers at several different altitudes. No well formulated treatment of this 
problem has yet been presented in the literature although techniques are cur- 
rently being developed by Smith (private communication). 


CONSTITUENT INVERSION 

Once we have obtained the temperature profile from an absorption band of 
a uniformly mixed atmospheric constituent, we can in principle infer the vertical 
distribution of nonuniformly mixed optically active gases which are in local 
thermodynamic equilibrium with the other atmospheric constituents. 

For purposes of formulating this problem, it is convenient to use the in- 
tegration by parts form of the transfer equation 

f X, dB 

I(^) = B(v, 0) + ^ [v, T( x)j r[u, u(x)] dx (24) 

Jo 


14 



where u(x) is the mass of the optically active gas between level x and the top of 
the atmosphere. Given measurements of I(v) and T(x), (24) is a nonlinear in- 
tegral equation for the unknown function u(x). 

Since r is generally a strong function of v unlike the Planck function in the 
temperature inversion problem, a somewhat different approach is required. 

One possible approach suggested by King (1964) is to do the temperature in- 
version in both a C0 2 band and an absorption band of the gas whose vertical 
distribution is to be determined. The former will give temperature as a func- 
tion of pressure and the latter will give temperature as a function of the ab- 
sorber mass of the gas whose distribution is unknown. The temperature can 
then be eliminated from the two relations to obtain absorber mass as a function 
of pressure. The possibility of applying this technique to measurements in the 
rotational bands of water vapor is currently being investigated by Yamamoto 
(private communication). Another method involving the iterative comparison 
of calculated and measured intensities in relatively broad spectral bands is 
being persued at ESSA (Smith, private communication). The object of the 
method is to obtain relatively crude two-parameter water vapor and tempera- 
ture distributions in the lower troposphere. 

Still another approach is being persued in our laboratory. The method in- 
volves an attempt at direct linearization of (24). It is assumed that the trans- 
mission function appearing in (24) can be approximated by 

r[v, u(x)l ^ t[v,u o (x)] + (fj) Su(x) (25) 

\ /u 0 (x) 


where u 0 (x) is some initial guess at the absorber distribution. When (25) is 
substituted into (24) we obtain an integral equation of the first kind in the 
correction function Su(x) 


I(*) - I„ ( v > 



K Q (v, x) 8u(x) dx 


(26) 


where I 0 {v) is the calculated intensity from the distribution u 0 (x) and the 
kernel is given by 

<?B (dr\ 

K o ( v - x ) ~ (*• , , 

\ /u 0 (x) 


(27) 
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The techniques developed for the temperature inversion problem can now be 
applied to (26). The process can be iterated as many times as desired. How- 
ever, when I(i') - I 0 (v) becomes smaller than the probable error of the meas- 
urements nothing further can be gained by iteration. 

Application of this method to the problem of obtaining water vapor distribu- 
tions from the 6.3m H 2 0 band has been made. The results of application of this 
method to synthetic data from a model atmosphere are shown in Figure 10. For 
application to the water vapor problem, 8 (log u*) is used rather than Su where 
u* is the reduced absorber mass. The curvature of the transmissivity as a 
function of log u* is small in the region of interest so even though the correction 
8 (log u*) may be relatively large, (25) is still a fairly good approximation. 

For this calculation, the transmissivities of Moller and Raschke (1964) were 
used. The temperature profile used in the model and the resulting kernel func- 
tions are shown in Figure 11. A two parameter water vapor distribution was 
assumed. In terms of the mixing ratio 



where q s and k are the parameters to be determined. 

The method has also been applied to data obtained with an IRIS instrument 
flown on a balloon at Palestine, Texas, on 8 May 1966. An average of thirteen 
spectra distributed throughout the day was employed. The inferred distribution, 
which can be regarded as an average over the day, is shown in Figure 12. 
Radiosonde measurements taken at the stations nearest Palestine are shown for 
comparison. 

Although only water vapor has been mentioned in detail here, the method 
should also be applicable to 0 3 in the earth's atmosphere. 

More sophisticated representations of u(x) than those indicated above can 
of course be employed. However, the more degrees of freedom u(x) is allowed, 
the more severe the stability problem will become. One obvious approach 
would be the application of empirical orthogonal functions. Whether or not 
water vapor is statistically stable enough to make the use of such functions 
practical is a subject of investigation. 
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SUMMARY 


The inverse problem of radiative transfer for planetary atmospheres has been 
reviewed with emphasis placed on the temperature and water vapor inversion 
problems for the terrestrial atmosphere. The general instability of the problem 
is the main factor to be dealt with in formulating a method of solution. No en- 
tirely satisfactory method has been developed for the solution of the problem. 
However, a few methods showing some promise have been formulated. 

The linear temperature inversion method developed by Wark and his as- 
sociates seems most promising for applications to the problem of the acquisition 
of vertical temperature profiles in the terrestrial atmosphere on a global basis. 
When used in conjunction with empirical orthogonal functions, the method allows 
the use of a maximum of a priori information on the vertical atmospheric 
structure. However, since the method is essentially statistical, we must expect 
results of a statistical nature. We cannot hope to recover profiles which may 
be strongly atypical. The success of the method will depend largely on our 
ability to estimate satisfactorily the temperature covariance. Best results can 
be expected for geographical areas and seasons for which the atmospheric 
structure is stable in a statistical sense. It would appear that the problem of 
calculating covariance functions requires further work. 

The only constraint imposed in the nonlinear method of King is that the 
temperature profile be represented by a series of ramps. It does not utilize 
additional a priori information, and the smoothing procedure is objective. Be- 
cause of the nonlinear nature of the method, error analysis is somewhat dif- 
ficult. The best approach seems to be by way of application of the method to a 
wide variety of climatological profiles. Such an analysis is currently being 
undertaken by King and his associates. Preliminary results indicate two ramp 
inferences can generally be expected. It would appear that this method may be 
most applicable to other planetary atmospheres for which a minimum amount 
of a priori information is available and only gross structure is sought. 

Recently the problem of inferring vertical water vapor distribution has 
been given some consideration. Work completed thus far indicates at least 
the gross features of the vertical distribution can be obtained from data taken 
with instruments such as the IRIS. The maximum amount of information which 
can be extracted is a subject of current investigations. No detailed analysis of 
the problem of ozone inversion has been made to our knowledge. However, 
methods developed for the water vapor inversion problem should be applicable, 
at least in principle. 
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The usefulness of information obtained on the vertical structure of the 
earth's atmosphere on a global basis will ultimately have to be determined by 
the meteorologists who wish to make use of the data. While analyses based 
on model atmospheres and data obtained from balloon flights at one or a few 
locations is useful in the development of techniques, only after satellite data 
covering a wide variety of seasons and geographic locations becomes available 
can the practicality of the approach be established. 
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• FOLDING PRESSURE VS WAVENUMBER 
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WAVENUMBER (crrr 1 ) 

Figure 1. Fifteen micron carbon dioxide absorption band e-folding pressure levels 
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Hgure 3. Temperature inversion of IRIS balloon flight data. The effects of the Lagrangian 

multiplier yon the inversion is demonstrated. 
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Figure 4. Two members of a set of emperical orthogonal functions. The function set 
was constructed from an ensemble of twenty-five temperature profiles 
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Figure 5. Synthetic data inversion in terms of empirical orthogonal functions. The 
sounding employed was a member of the ensemble used in constructing the function set 
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Figure 6. Synthetic data inversion in terms of empirical orthogonal functions. The sounding 
employed was not a member of the ensemble used in constructing the function set 
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NONLINEAR INVERSION OF SYNTHETIC DATA 
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Figure 7. An application of the nonlinear inversion method to synthetic data. This illustrates the 
behavior of the method when applied to a tropical profile 
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thod to synthetic data. The behavior of the 
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Figure 9. Nonlinear inversion of a model Martian atmosphere. IRIS-type data was 

synthesized numerically 
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RELATIVE WEIGHT TEMPERATURE (°K) 

Figure 11. Weighting functions employed in the water vapor inversion shown in Figure 9. The temperature 

profile employed in the model is also shown 
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MIXING RATIO (G/KG) 


